# load necessary packages and data
library(tidyr)
library(dplyr)
library(readr)
library(plotly)
library(ggplot2)
library(jgcricolors)
library(ggsci)
library(Polychrome)

# set directories
run_dir <- "C:/Users/spei632/Documents/GCAM_industry/food_processing/initial_breakout_results/data/food_processing_validation_run_11-9-23"
fig_dir <- paste0(run_dir,"/figures_validation")
if (!dir.exists(fig_dir)) {dir.create(fig_dir)}
results_dir <- fig_dir
setwd(run_dir)

# load functions
source("C:/Users/spei632/Documents/GCAM_industry/food_processing/initial_breakout_results/parse_PIC_data.R")

# load data
queries <- list.files(run_dir, pattern = 'queryoutall')
for (i in queries) {
  filename <- gsub('.csv','', i) %>% gsub('queryoutall_','', .)
  if (!grepl("food_pro", filename) & !grepl("national_account", filename)) {
  # if (!grepl("food_pro", filename)) {  
    assign(filename, read_csv(paste0(run_dir, "/", i), skip = 1) %>% 
             parse_PIC_data() %>%
             # indicate if master or food processing branch
             mutate(branch = if_else(grepl("date=2023-9-11", date), "food_industry", "master"),
                    scenario = if_else(scenario == "Reference", "GCAM", scenario),
                    policy_type = case_when(grepl("2p6", scenario) ~"RCP 2p6", 
                                            grepl("Tax25", scenario) ~ "Tax25",
                                            TRUE ~ "No policy"),
                    scenario_type = case_when(grepl("SSP", scenario) ~ substr(scenario, 6, 9),
                                              grepl("USA", scenario) ~ scenario,
                                              TRUE ~ "GCAM")))
  } else {
    # if it is a food processing query, there will be no results for some of the scenarios, and we need to handle separately
    # figure out which line to skip to - first read in normally
    temp <- read_csv(paste0(run_dir, "/", i), col_names = FALSE)
    # get first index of when actual data starts
    index_start <- which(grepl("scenario", temp$X1))[1]
    # get the comment string indicating no values were output, to skip this in reading csv
    comment_string <- temp$X1[which(grepl("had error:", temp$X1))[1]]
    assign(filename, read_csv(paste0(run_dir, "/", i), skip = index_start - 1, comment = comment_string) %>% 
             parse_PIC_data() %>%
             # indicate if master or food processing branch
             mutate(branch = if_else(grepl("date=2023-9-11", date), "food_industry", "master"),
                    scenario = if_else(scenario == "Reference", "GCAM", scenario),
                    policy_type = case_when(grepl("2p6", scenario) ~"RCP 2p6", 
                                            grepl("Tax25", scenario) ~ "Tax25",
                                            TRUE ~ "No policy"),
                    scenario_type = case_when(grepl("SSP", scenario) ~ substr(scenario, 6, 9),
                                              grepl("USA", scenario) ~ scenario,
                                              TRUE ~ "GCAM")))
    rm(temp)
    rm(index_start)
    rm(comment_string)
  }
}

# load GTAP data
CostShare_FoodProcessing_GTAP <- read_csv("C:/Users/spei632/Documents/GCAM_industry/food_processing/initial_exploration/data/CostShare_FoodProcessing_GTAP.csv")

# mappings
GWP_AR5 <- read_csv("C:/Users/spei632/Documents/GCAM_industry/food_processing/mappings/GWP_AR5.csv")
industry_sector_agg <- read_csv("C:/Users/spei632/Documents/GCAM_industry/food_processing/mappings/industry_sector_agg.csv")

# constants
C_to_CO2 <- 44/12

# get color palettes
pal_sel <- jgcricol()$pal_all
pal_sel_extended <- c(pal_sel, "food processing" = "mediumpurple", "biomass cogen" = "seagreen3",
                      "coal cogen" = "dimgray", "gas cogen" = "skyblue2", "gas with solar" = "lightblue",
                      "refined liquids cogen" = "lightpink", "district heat" = "darkslategray",
                      "process heat food processing" = "mediumpurple1", "electricity with solar" = "antiquewhite2", "electric heat pump" = "lavenderblush3")
data(palette36)
palette36 <- unname(palette36)
pal_industry_sectors <- c("food processing" = "red", "agricultural energy use" = "burlywood3", "aluminum" = "lightsalmon", "cement" = "khaki2", "chemical" = "lightblue", "construction" = "plum", "iron and steel" = "palegreen3", "mining energy use" = "tan3", "N fertilizer" = "pink", "other industry" = "lightgray", "paper" = "wheat")
pal_scenarios <- c("GCAM" = "red", "GCAM_2p6" = "lightpink","GCAM_SSP1" = "tan3", "GCAM_SSP1_2p6" = "burlywood3", "GCAM_SSP2" = "mediumpurple1", "GCAM_SSP2_2p6" = "plum", "GCAM_SSP3" = "khaki2", "GCAM_SSP4" = "skyblue2", "GCAM_SSP4_2p6" = "lightblue", "GCAM_SSP5" = "palegreen3", "GCAM_SSP5_2p6" = "darkseagreen")
pal_scenarios_w_USA <- c(pal_scenarios, "GCAM-USA_Ref" = "lavenderblush3", "GCAM-USA_Tax25" = "lightgray")

# scenarios to select
ref_scenario_name <- "GCAM"
food_branch_name <- "food_industry"
master_branch_name <- "master"

# industry sectors
industry_sectors <- unique(industry_shrwts_subsector$sector)

Food processing sector breakout branch results

Reference scenario

Food processing sector outputs

Food processing output - just confirming that this matches food demand.

# aggregate food demand to total calories
food_demand_total <- food_demand %>%
  group_by(scenario, branch, policy_type, scenario_type, Units, region, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

food_demand_total_food_pro_total_output <- food_demand_total %>%
  mutate(sector = "food demand") %>%
  bind_rows(food_pro_prod_region %>% dplyr::select(Units, scenario, branch, policy_type, scenario_type, region, year, sector, value))

ggplot(food_demand_total_food_pro_total_output %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, color = region)) +
  geom_line() +
  facet_wrap(~sector, ncol = 2, scale = "free_y") +
  scale_color_manual(values = palette36) +
  labs(x = "", y = "Pcal", title = "Food processing and food demand output by region, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank(),
        legend.position = "bottom")

ggsave(paste0(fig_dir, "/food_pro_food_demand_ref_food_ind.png"), width = 20, height = 15)

Process heat food processing output by technology.

food_pro_process_heat_prod_tech_heat_only <- food_pro_process_heat_prod_tech %>%
  filter(output == "process heat food processing")

# get total
food_pro_process_heat_prod_heat_only_total <- food_pro_process_heat_prod_tech_heat_only %>%
  group_by(scenario, region, sector, Units, year, branch, policy_type, scenario_type) %>%
  summarize(value = sum(value)) %>%
  ungroup()

# get total global
food_pro_process_heat_prod_heat_only_total_global <- food_pro_process_heat_prod_heat_only_total %>%
  group_by(scenario, branch, policy_type, scenario_type, sector, Units, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

ggplot(food_pro_process_heat_prod_tech_heat_only %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = technology)) +
  geom_col() +
  facet_wrap(~region, ncol = 8, scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Outputs of process heat food processing by technology, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_by_tech_ref_food_ind.png"), width = 30, height = 20)

ggplot(food_pro_process_heat_prod_tech_heat_only %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = technology)) +
  geom_col() +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Global outputs of process heat food processing by technology\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_by_tech_ref_food_ind_global.png"), width = 10, height = 8)

Food processing sector energy inputs

Total food processing energy use.

food_pro_energy_inputs_edit <- food_pro_inputs %>%
  filter(!grepl("water", input)) %>%
  mutate(input_edit = case_when(input == "elect_td_ind" ~ "electricity",
                                input == "refined liquids industrial" ~ "refined liquids",
                                input == "delivered biomass" ~ "biomass",
                                input == "delivered coal" ~ "coal",
                                input == "wholesale gas" ~ "gas",
                                input == "global solar resource" ~ "solar",
                                TRUE ~ input))

ggplot(food_pro_energy_inputs_edit %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_wrap(~region, ncol = 8, scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Energy inputs to food processing and process heat food processing, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_ref_food_ind.png"), width = 30, height = 20)

ggplot(food_pro_energy_inputs_edit %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Global energy inputs to food processing and process heat food processing\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_ref_food_ind_global.png"), width = 10, height = 8)

# separate by process heat food processing and food processing
food_pro_inputs_process_heat_sep_edit <- food_pro_inputs_process_heat_sep %>%
  filter(!grepl("water", input)) %>%
  mutate(input_edit = case_when(input == "elect_td_ind" ~ "electricity",
                                input == "refined liquids industrial" ~ "refined liquids",
                                input == "delivered biomass" ~ "biomass",
                                input == "delivered coal" ~ "coal",
                                input == "wholesale gas" ~ "gas",
                                input == "global solar resource" ~ "solar",
                                TRUE ~ input))

ggplot(food_pro_inputs_process_heat_sep_edit %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  facet_wrap(~sector, ncol = 2) +
  labs(x = "", y = "EJ", title = "Global energy inputs to food processing and process heat food processing\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_sep_sector_ref_food_ind_global.png"), width = 15, height = 8)

# plot just for 8 largest energy users in 2100
top_8_energy_users_food_pro_2100_ref <- food_pro_energy_inputs_edit %>%
  filter(year == 2100 & scenario == "GCAM") %>%
  group_by(region) %>%
  summarize(value = sum(value)) %>%
  arrange(desc(value)) %>%
  pull(region) %>%
  head(8)
  
ggplot(food_pro_energy_inputs_edit %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005 & region %in% top_8_energy_users_food_pro_2100_ref) %>%
         mutate(region = factor(region, levels = top_8_energy_users_food_pro_2100_ref)),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_wrap(~region, nrow = 4) +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Food processing energy inputs, reference scenario", 
       caption = "Top 8 food processing energy consumers in 2100 are shown") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_ref_food_ind_top_8_en_users_2100.png"), width = 8, height = 10)

Food processing overall sector inputs.

food_pro_energy_inputs_overall_sector <- inputs_tech %>%
  filter(sector == "food processing" & !grepl("water", input)) %>%
  mutate(input = if_else(input == "elect_td_ind", "electricity", input))

ggplot(food_pro_energy_inputs_overall_sector %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input)) +
  geom_col() +
  facet_wrap(~region, ncol = 8, scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Energy inputs to food processing, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_overall_sector_energy_inputs_ref_food_ind.png"), width = 30, height = 20)

ggplot(food_pro_energy_inputs_overall_sector %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input)) +
  geom_col() +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Global energy inputs to food processing\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_overall_sector_energy_inputs_ref_food_ind_global.png"), width = 10, height = 8)

Process heat food processing energy inputs.

food_pro_process_heat_energy_inputs_edit <- food_pro_inputs_process_heat_sep_edit %>%
  filter(sector == "process heat food processing")

ggplot(food_pro_process_heat_energy_inputs_edit %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_wrap(~region, ncol = 8, scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Energy inputs to process heat food processing, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_energy_inputs_ref_food_ind.png"), width = 30, height = 20)

ggplot(food_pro_process_heat_energy_inputs_edit %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Global energy inputs to process heat food processing\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_energy_inputs_ref_food_ind_global.png"), width = 10, height = 8)

Food processing energy use coefficient per calorie consumed. The coefficient based on electricity and process heat is calculated as the sum of electricity + process heat input into overall food processing, while the coefficient based on input energy is calculated as the sum of the energy inputs into process heat food processing and overall food processing.

# get total energy use
food_pro_energy_total <- food_pro_energy_inputs_edit %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year) %>%
  summarize(EJ = sum(value)) %>%
  ungroup()

# get electricity + process heat inputs
food_pro_energy_elec_process_total <- food_pro_energy_inputs_overall_sector %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year) %>%
  summarize(EJ_elec_process = sum(value)) %>%
  ungroup()

# get coefficient of energy use per calorie consumed 
food_pro_coef <- food_pro_energy_total %>%
  left_join(food_pro_energy_elec_process_total) %>%
  left_join(food_demand_total %>%
              dplyr::select(scenario, branch, policy_type, scenario_type, region, year, Pcal = value)) %>%
  mutate(coef_input_energy = EJ / Pcal,
         coef_electricity_process_heat = EJ_elec_process / Pcal) %>%
  dplyr::select(-c(EJ, EJ_elec_process, Pcal)) %>%
  pivot_longer(c(coef_input_energy, coef_electricity_process_heat), names_to = "coef_type", values_to = "value")

ggplot(food_pro_coef %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, color = coef_type)) +
  geom_line(size = 1) + 
  facet_wrap(~region, ncol = 8) +
  labs(x = "", y = "EJ per Pcal", title = "Food processing energy use coefficient, reference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_energy_use_coef_ref_food_ind.png"), width = 30, height = 20)

Food processing sector carbon emissions

# convert all emissions to CO2 from C
CO2_emis_sector <- CO2_emis_sector %>%
  mutate(Mt_CO2 = value * C_to_CO2)
CO2_emis_sector_no_bio <- CO2_emis_sector_no_bio %>% 
  mutate(Mt_CO2 = value * C_to_CO2)
CO2_emis_region <- CO2_emis_region %>%
  mutate(Mt_CO2 = value * C_to_CO2)
LUC_emis_region <- LUC_emis_region  %>%
  mutate(Mt_CO2 = value * C_to_CO2)

# get global
CO2_emis_sector_global <- CO2_emis_sector %>%
  group_by(scenario, branch, policy_type, scenario_type, sector, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()
CO2_emis_sector_no_bio_global <- CO2_emis_sector_no_bio %>%
  group_by(scenario, branch, policy_type, scenario_type, sector, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()
CO2_emis_global <- CO2_emis_region %>%
  group_by(scenario, branch, policy_type, scenario_type, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()
LUC_emis_global <- LUC_emis_region %>%
  group_by(scenario, branch, policy_type, scenario_type, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()

# get total CO2 emissions
CO2_emis_global_w_LUC <- CO2_emis_global %>%
  mutate(emis_type = "FFI") %>%
  rbind(LUC_emis_global %>% 
          mutate(emis_type = "LUC") %>%
          filter(year %in% unique(CO2_emis_global$year)))

CO2_emis_global_w_LUC_total <- CO2_emis_global_w_LUC %>%
  group_by(scenario, branch, policy_type, scenario_type, year) %>% 
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()

# note that because the food processing overall sector takes in only process heat food processing and electricity, it has no associated direct emissions
# all food processing emissions come from the process heat food processing sector
ggplot(CO2_emis_sector_no_bio %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005 & sector == "process heat food processing"),
       aes(x = year, y = Mt_CO2)) +
  geom_line(size = 1) + 
  facet_wrap(~region, nrow = 4) + 
  labs(x = "", y = "Mt CO2", title = "Food processing CO2 emissions (no bio), reference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_food_pro_ref_food_ind.png"), width = 30, height = 20)

ggplot(CO2_emis_sector_no_bio_global %>% 
         filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005 & sector == "process heat food processing"),
       aes(x = year, y = Mt_CO2)) +
  geom_line(size = 1) + 
  labs(x = "", y = "Mt CO2", title = "Global food processing CO2 emissions (no bio)\nReference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_food_pro_ref_food_ind_global.png"), width = 10, height = 8)

Industry energy use by sector

# aggregate industry sectors
industry_final_en_sector_agg <- industry_final_en_tech_fuel %>%
  filter(sector %in% industry_sectors) %>%
  left_join(industry_sector_agg) %>%
  group_by(Units, scenario, branch, policy_type, scenario_type, region, sector_agg, year) %>%
  summarize(value = sum(value)) %>%
  ungroup() %>%
  mutate(sector_agg = factor(sector_agg, 
                             levels = c("food processing", unique((industry_sector_agg %>% filter(sector_agg != "food processing"))$sector_agg))))

industry_final_en_sector_agg_global <- industry_final_en_sector_agg %>%
  group_by(Units, scenario, branch, policy_type, scenario_type, sector_agg, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

# aggregate to total industry
industry_final_en_sector_agg_total <- industry_final_en_sector_agg %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year, Units) %>%
  summarize(value = sum(value)) %>%
  ungroup()

industry_final_en_sector_agg_total_global <- industry_final_en_sector_agg_total %>%
  group_by(scenario, branch, policy_type, scenario_type, year, Units) %>%
  summarize(value = sum(value)) %>%
  ungroup()

ggplot(industry_final_en_sector_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = sector_agg)) +
  geom_col() +
  facet_wrap(~region, nrow = 4) +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Industry final energy use by sector, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_sector_ref_food_ind.png"), width = 30, height = 20)

ggplot(industry_final_en_sector_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = sector_agg)) +
  geom_col() +
  facet_wrap(~region, nrow = 4, scales = "free_y") +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Industry final energy use by sector, reference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_sector_ref_food_ind_free_y.png"), width = 30, height = 20)

ggplot(industry_final_en_sector_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = value, fill = sector_agg)) +
  geom_col() +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Global industry final energy use by sector\nReference scenario, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_sector_ref_food_ind_global.png"), width = 10, height = 8)

Industry carbon emissions by sector

# aggregate sectors
CO2_emis_sector_no_bio_industry_agg <- CO2_emis_sector_no_bio %>%
  filter(sector %in% industry_sectors) %>%
  left_join(industry_sector_agg) %>%
  group_by(scenario, branch, policy_type, scenario_type, region, sector_agg, year) %>%
  summarize(value = sum(value), Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup() %>%
  mutate(sector_agg = factor(sector_agg, 
                             levels = c("food processing", unique((industry_sector_agg %>% filter(sector_agg != "food processing"))$sector_agg))))

CO2_emis_sector_no_bio_industry_agg_global <- CO2_emis_sector_no_bio_industry_agg  %>%
  group_by(scenario, branch, policy_type, scenario_type, sector_agg, year) %>%
  summarize(value = sum(value),
            Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()

# aggregate to total industry
CO2_emis_sector_no_bio_industry_agg_total <- CO2_emis_sector_no_bio_industry_agg %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()

CO2_emis_sector_no_bio_industry_agg_total_global <- CO2_emis_sector_no_bio_industry_agg_total %>%
  group_by(scenario, branch, policy_type, scenario_type, year) %>%
  summarize(Mt_CO2 = sum(Mt_CO2)) %>%
  ungroup()
  
ggplot(CO2_emis_sector_no_bio_industry_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = Mt_CO2, fill = sector_agg)) +
  geom_col() + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2", title = "Industry CO2 emissions (no bio), reference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_sector_agg_ref_food_ind.png"), width = 30, height = 20)

ggplot(CO2_emis_sector_no_bio_industry_agg_global %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = Mt_CO2, fill = sector_agg)) +
  geom_col() + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  labs(x = "", y = "Mt CO2", title = "Global industry CO2 emissions (no bio)\nReference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_sector_agg_ref_food_ind_global.png"), width = 10, height = 8)

Industry non-CO2 emissions by sector

# get industry non-CO2 emissions
nonCO2_emis_sector_CO2eq_industry <- nonCO2_emis_sector %>%
  filter(sector %in% industry_sectors) %>%
  left_join(GWP_AR5 %>% rename(GHG = ghg)) %>%
  filter(!is.na(GWP)) %>%
  mutate(Mt_CO2e = value * GWP)

nonCO2_emis_sector_total_CO2eq_industry <- nonCO2_emis_sector_CO2eq_industry %>%
  group_by(scenario, branch, policy_type, scenario_type, region, sector, year) %>%
  summarize(Mt_CO2e = sum(Mt_CO2e)) %>%
  ungroup()

# aggregate sectors
nonCO2_emis_sector_total_CO2eq_industry_agg <- nonCO2_emis_sector_total_CO2eq_industry %>%
  left_join(industry_sector_agg) %>%
  group_by(scenario, branch, policy_type, scenario_type, region, sector_agg, year) %>%
  summarize(Mt_CO2e = sum(Mt_CO2e)) %>%
  ungroup() %>%
  mutate(sector_agg = factor(sector_agg, 
                             levels = c("food processing", unique((industry_sector_agg %>% filter(sector_agg != "food processing"))$sector_agg))))

# aggregate to total industry
nonCO2_emis_sector_total_CO2eq_industry_agg_total <- nonCO2_emis_sector_total_CO2eq_industry_agg %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year) %>%
  summarize(Mt_CO2e = sum(Mt_CO2e)) %>%
  ungroup()

nonCO2_emis_sector_total_CO2eq_industry_agg_total_global <- nonCO2_emis_sector_total_CO2eq_industry_agg_total %>%
  group_by(scenario, branch, policy_type, scenario_type, year) %>%
  summarize(Mt_CO2e = sum(Mt_CO2e)) %>%
  ungroup()

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = Mt_CO2e, fill = sector_agg)) +
  geom_col() + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2e", title = "Industry GHG emissions, reference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_sector_ref_food_ind.png"), width = 30, height = 20)

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005),
       aes(x = year, y = Mt_CO2e, fill = sector_agg)) +
  geom_col() + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  labs(x = "", y = "Mt CO2e", title = "Global industry GHG emissions\nReference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_sector_ref_food_ind_global.png"), width = 10, height = 8)

Industry water withdrawals by sector

# water withdrawals
industry_water_withdrawals <- industry_water_withdrawals_consumption %>% 
  filter(input == "water_td_ind_W")

industry_water_withdrawals_global <- industry_water_withdrawals %>%
  filter(region != "Global") %>%
  group_by(Units, scenario, branch, policy_type, scenario_type, sector, input, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

# aggregate to total excluding H2 sectors
industry_water_withdrawals_total_no_H2 <- industry_water_withdrawals %>%
  filter(!grepl("H2", sector)) %>%
  group_by(scenario, branch, policy_type, scenario_type, region, year, Units) %>%
  summarize(value = sum(value)) %>%
  ungroup()

industry_water_withdrawals_total_no_H2_global <- industry_water_withdrawals_total_no_H2 %>%
  filter(region != "Global") %>%
  group_by(scenario, branch, policy_type, scenario_type, year, Units) %>%
  summarize(value = sum(value)) %>%
  ungroup()

ggplot(industry_water_withdrawals %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005 & !grepl("H2", sector)),
       aes(x = year, y = value, fill = sector)) +
  geom_col() +
  facet_wrap(~region, ncol = 4, scales = "free_y") + 
  scale_fill_manual(values = pal_industry_sectors, limits = force, drop = TRUE) + 
  labs(x = "", y = "km^3", title = "Industry water withdrawals, reference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_ref_food_ind.png"), width = 30, height = 20)

ggplot(industry_water_withdrawals_global %>% filter(scenario == ref_scenario_name & branch == food_branch_name & year >= 2005 & !grepl("H2", sector)),
       aes(x = year, y = value, fill = sector)) +
  geom_col() +
  scale_fill_manual(values = pal_industry_sectors, limits = force, drop = TRUE) + 
  labs(x = "", y = "km^3", title = "Global industry water withdrawals\nReference scenario, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_ref_food_ind_global.png"), width = 10, height = 8)

SSP/RCP results

Food processing total energy use

ggplot(industry_final_en_sector_agg %>% filter(branch == food_branch_name & year >= 2005 & sector_agg == "food processing" & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "EJ", title = "Food processing total energy use, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_final_en_scenario_comp.png"), width = 30, height = 20)

ggplot(industry_final_en_sector_agg_global %>% filter(branch == food_branch_name & year >= 2005 & sector_agg == "food processing" & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "EJ", title = "Global food processing total energy use, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_final_en_scenario_comp_global.png"), width = 15, height = 8)

Food processing energy use by fuel

food_pro_energy_inputs_edit_global <- food_pro_energy_inputs_edit %>%
  group_by(scenario, branch, policy_type, scenario_type, sector, input_edit, Units, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

ggplot(food_pro_energy_inputs_edit %>% filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_grid(rows = vars(policy_type), cols = vars(scenario_type), scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Global food processing total energy use by fuel, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_scenario_comp_global.png"), width = 18, height = 12)
  
ggplot(food_pro_energy_inputs_edit_global %>% filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) +
  scale_color_manual(values = pal_scenarios) + 
  facet_grid(rows = vars(policy_type), cols = vars(input_edit), scale = "free_y") +
  labs(x = "", y = "EJ", title = "Global food processing total energy use by fuel, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_by_fuel_scenario_comp_global.png"), width = 18, height = 10)

# plot just for 8 largest energy users in 2100 in 2.6 scenario
top_8_energy_users_food_pro_2100_2p6 <- food_pro_energy_inputs_edit %>%
  filter(year == 2100 & scenario == "GCAM_2p6") %>%
  group_by(region) %>%
  summarize(value = sum(value)) %>%
  arrange(desc(value)) %>%
  pull(region) %>%
  head(8)
  
ggplot(food_pro_energy_inputs_edit %>% filter(scenario == "GCAM_2p6" & branch == food_branch_name & year >= 2005 & region %in% top_8_energy_users_food_pro_2100_2p6) %>%
         mutate(region = factor(region, levels = top_8_energy_users_food_pro_2100_ref)),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_wrap(~region, nrow = 4) +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Food processing energy inputs, RCP 2.6 scenario", 
       caption = "Top 8 food processing energy consumers in 2100 are shown") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_total_energy_inputs_2p6_food_ind_top_8_en_users_2100.png"), width = 8, height = 10)

Process heat food processing outputs by technology

ggplot(food_pro_process_heat_prod_tech_heat_only %>% 
         filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, fill = technology)) +
  geom_col() +
  facet_grid(rows = vars(policy_type), cols = vars(scenario_type), scale = "free_y") +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Outputs of process heat food processing by technology, food industry branch") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_by_tech_scenario_comp_global.png"), width = 18, height = 12)

ggplot(food_pro_process_heat_prod_heat_only_total %>% filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "EJ", title = "Process heat food processing production, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_total_scenario_comp.png"), width = 30, height = 20)

ggplot(food_pro_process_heat_prod_heat_only_total_global %>% filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "EJ", title = "Global process heat food processing production, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_total_scenario_comp_global.png"), width = 15, height = 8)

# plot just the main GCAM scenarios
ggplot(food_pro_process_heat_prod_tech_heat_only %>% 
         filter(branch == food_branch_name & year >= 2005 & scenario %in% c("GCAM", "GCAM_2p6")) %>%
         mutate(name = if_else(scenario == "GCAM", "Reference", "RCP 2.6"),
                name = factor(name, levels = c("Reference", "RCP 2.6"))),
       aes(x = year, y = value, fill = technology)) +
  geom_col() +
  facet_wrap(~name) +
  scale_fill_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "technology") +
  labs(x = "", y = "EJ", title = "Global food processing process heat production by technology") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_process_heat_outputs_by_tech_scenario_comp_global_sel.png"), width = 9, height = 6)

Food processing carbon emissions

ggplot(CO2_emis_sector_no_bio %>% 
         filter(branch == food_branch_name & year >= 2005 & sector == "process heat food processing" & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4) + 
  labs(x = "", y = "Mt CO2", title = "Food processing CO2 emissions (no bio), food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_food_pro_scenario_comp.png"), width = 30, height = 20)

ggplot(CO2_emis_sector_no_bio_global %>% 
         filter(branch == food_branch_name & year >= 2005 & sector == "process heat food processing" & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) +
  labs(x = "", y = "Mt CO2", title = "Global food processing CO2 emissions (no bio), food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_food_pro_scenario_comp_global.png"), width = 10, height = 8)

Food processing prices

ggplot(food_pro_prices %>% filter(branch == food_branch_name & year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, color = scenario)) +
  geom_line(size = 1) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4) + 
  labs(x = "", y = "1975$/Mcal", title = "Food processing prices, food industry branch") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/food_pro_prices_scenario_comp.png"), width = 30, height = 20)

Comparison with master branch

Reference scenario

Industry energy use

Total industry energy use.

# compare total industry energy use across regions and globally
ggplot(industry_final_en_sector_agg_total %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "EJ", title = "Total industry energy use, reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/total_industry_energy_use_ref_comp.png"), width = 30, height = 20)

ggplot(industry_final_en_sector_agg_total_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "EJ", title = "Global total industry energy use, reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/total_industry_energy_use_ref_comp_global.png"), width = 10, height = 8)

Industry energy use by sector.

fig <- ggplot(industry_final_en_sector_agg %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, fill = sector_agg)) + 
  geom_col() +
  facet_grid(rows = vars(region), cols = vars(branch), scales = "free_y") +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Industry final energy use by sector, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())
ggsave(paste0(fig_dir, "/industry_final_en_sector_ref_comp.png"), fig, width = 10, height = 30)

ggplot(industry_final_en_sector_agg %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, fill = sector_agg)) + 
  geom_col() +
  facet_wrap(~branch, ncol = 2) +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Global industry final energy use by sector, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_sector_ref_comp_global.png"), width = 15, height = 8)

Industry energy use by fuel.

industry_final_en_fuel <- industry_final_en_tech_fuel %>%
  group_by(scenario, branch, policy_type, scenario_type, region, input, Units, year) %>%
  summarize(value = sum(value)) %>%
  ungroup() %>%
  mutate(input_edit = case_when(input == "elect_td_ind" ~ "electricity",
                                input == "refined liquids industrial" ~ "refined liquids",
                                input == "delivered biomass" ~ "biomass",
                                input == "delivered coal" ~ "coal",
                                input == "wholesale gas" ~ "gas",
                                input == "global solar resource" ~ "solar",
                                grepl("H2", input) ~ "hydrogen",
                                TRUE ~ input))

industry_final_en_fuel_global <- industry_final_en_fuel %>%
  group_by(scenario, branch, policy_type, scenario_type, input, input_edit, Units, year) %>%
  summarize(value = sum(value)) %>%
  ungroup()

ggplot(industry_final_en_fuel %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, fill = input_edit)) +
  geom_col() +
  facet_wrap(~branch, ncol = 2, scale = "free_y") +
  scale_fill_manual(values = c(pal_sel_extended, "regional woodpulp for energy" = "darkgreen"), drop = TRUE, limits = force, name = "fuel") +
  labs(x = "", y = "EJ", title = "Global industry final energy use by fuel, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_fuel_ref_comp_global.png"), width = 15, height = 8)

ggplot(industry_final_en_fuel_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch)) +
  geom_line(size = 1) +
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~input, ncol = 4, scale = "free_y") +
  labs(x = "", y = "EJ", title = "Global industry final energy use by fuel, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_fuel_by_fuel_ref_comp_global.png"), width = 15, height = 8)

Industry carbon emissions

ggplot(CO2_emis_sector_no_bio_industry_agg_total %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2", title = "Total industry CO2 emissions (no bio), reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_total_ref_comp.png"), width = 30, height = 20)

ggplot(CO2_emis_sector_no_bio_industry_agg_total_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "Mt CO2", title = "Global total industry CO2 emissions (no bio), reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_total_ref_comp_global.png"), width = 10, height = 8)

fig <- ggplot(CO2_emis_sector_no_bio_industry_agg %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, fill = sector_agg)) +
  geom_col() + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  facet_grid(rows = vars(region), cols = vars(branch), scales = "free_y") + 
  labs(x = "", y = "Mt CO2", title = "Industry CO2 emissions (no bio), reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())
ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_sector_agg_region_ref_comp.png"), fig, width = 10, height = 30)

ggplot(CO2_emis_sector_no_bio_industry_agg_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, fill = sector_agg)) +
  geom_col() +
  facet_wrap(~branch, ncol = 2) + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  labs(x = "", y = "Mt CO2", title = "Global industry CO2 emissions (no bio), reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_sector_agg_ref_comp_global.png"), width = 10, height = 8)

Industry non-CO2 emissions

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg_total %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2e, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2e", title = "Total industry GHG emissions, reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_total_ref_comp.png"), width = 30, height = 20)

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg_total_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2e, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "Mt CO2e", title = "Global total industry GHG emissions, reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_total_ref_comp_global.png"), width = 10, height = 8)

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2e, fill = sector_agg)) +
  geom_col() + 
  facet_wrap(~branch, ncol = 2) + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  labs(x = "", y = "Mt CO2e", title = "Global industry GHG emissions, reference scenario") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_sector_ref_comp_global.png"), width = 15, height = 8)

Industry water withdrawals

ggplot(industry_water_withdrawals_total_no_H2 %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "km^3", title = "Total industry water withdrawals, reference scenario\n(excluding hydrogen sectors)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_total_ref_comp.png"), width = 30, height = 20)

ggplot(industry_water_withdrawals_total_no_H2_global %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "km^3", title = "Global total industry water withdrawals, reference scenario\n(excluding hydrogen sectors)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_total_ref_comp_global.png"), width = 10, height = 8)

fig <- ggplot(industry_water_withdrawals %>% filter(scenario == ref_scenario_name & year >= 2005 & !grepl("H2", sector)),
       aes(x = year, y = value, fill = sector)) +
  geom_col() +
  facet_grid(rows = vars(region), cols = vars(branch), scales = "free_y") +
  scale_fill_manual(values = pal_industry_sectors, limits = force, drop = TRUE) +
  labs(x = "", y = "km^3", title = "Industry water withdrawals, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())
ggsave(paste0(fig_dir, "/industry_water_withdrawals_ref_comp.png"), fig, width = 10, height = 30)

ggplot(industry_water_withdrawals %>% filter(scenario == ref_scenario_name & year >= 2005 & !grepl("H2", sector)),
       aes(x = year, y = value, fill = sector)) +
  geom_col() +
  facet_wrap(~branch, ncol = 2) +
  scale_fill_manual(values = pal_industry_sectors, limits = force, drop = TRUE) +
  labs(x = "", y = "km^3", title = "Global industry water withdrawals, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_ref_comp_global.png"), width = 15, height = 8)

Industry fuel prices

industry_fuel_prices_edit <- industry_fuel_prices %>%
  mutate(fuel_edit = case_when(fuel == "elect_td_ind" ~ "electricity",
                               fuel == "refined liquids industrial" ~ "refined liquids",
                               fuel == "delivered biomass" ~ "biomass",
                               fuel == "delivered coal" ~ "coal",
                               fuel == "wholesale gas" ~ "gas",
                               grepl("H2", fuel) ~ "hydrogen",
                               TRUE ~ fuel))

ggplot(industry_fuel_prices_edit %>% filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = value, linetype = branch, color = fuel_edit)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_sel_extended, drop = TRUE, limits = force, name = "fuel") +
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "1975$/GJ", title = "Fuel prices to industry, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_fuel_prices_ref_comp.png"), width = 30, height = 20)

Total carbon emissions

ggplot(CO2_emis_global_w_LUC_total %>%  filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, linetype = branch)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "Mt CO2", title = "Global total net CO2 emissions, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_net_ref_comp_global.png"), width = 10, height = 8)

ggplot(CO2_emis_global %>%  filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, linetype = branch)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  labs(x = "", y = "Mt CO2", title = "Global total FFI CO2 emissions, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_FFI_ref_comp_global.png"), width = 10, height = 8)
  

ggplot(CO2_emis_region %>%  filter(scenario == ref_scenario_name & year >= 2005),
       aes(x = year, y = Mt_CO2, linetype = branch)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2", title = "Total FFI CO2 emissions, reference scenario") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_FFI_ref_comp.png"), width = 30, height = 20)

SSP/RCP results

Industry energy use

ggplot(industry_final_en_sector_agg_total %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "EJ", title = "Total industry energy use") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/total_industry_energy_use_scenario_branch_comp.png"), width = 30, height = 20)

ggplot(industry_final_en_sector_agg_total_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  facet_wrap(~policy_type, ncol = 2) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  labs(x = "", y = "EJ", title = "Global total industry energy use") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/total_industry_energy_use_scenario_branch_comp_global.png"), width = 15, height = 8)

ggplot(industry_final_en_fuel_global %>% 
         group_by(scenario, branch, policy_type, scenario_type, input_edit, Units, year) %>%
         summarize(value = sum(value)) %>%
         ungroup() %>%
         filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, linetype = branch, color = input_edit)) +
  geom_line(size = 1) +
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = c(pal_sel_extended, "regional woodpulp for energy" = "darkgreen"), drop = TRUE, limits = force, name = "fuel") + 
  facet_grid(rows = vars(policy_type), cols = vars(scenario_type), scale = "free_y") +
  labs(x = "", y = "EJ", title = "Global industry total energy use by fuel") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_fuel_scenario_branch_comp_global.png"), width = 18, height = 12)

ggplot(industry_final_en_sector_agg_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, fill = sector_agg)) +
  geom_col() +
  facet_grid(rows = vars(branch), cols = vars(scenario), scales = "free_y") +
  scale_fill_manual(values = pal_industry_sectors, name = "sector") +
  labs(x = "", y = "EJ", title = "Global industry total energy use by sector") +
  theme_bw() +
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_final_en_sector_scenario_branch_comp_global.png"), width = 30, height = 12)

Industry carbon emissions

ggplot(CO2_emis_sector_no_bio_industry_agg_total %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2", title = "Total industry CO2 emissions (no bio)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_total_scenario_branch_comp.png"), width = 30, height = 20)

ggplot(CO2_emis_sector_no_bio_industry_agg_total_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "Mt CO2", title = "Global total industry CO2 emissions (no bio)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_total_scenario_branch_comp_global.png"), width = 15, height = 8)

ggplot(CO2_emis_sector_no_bio_industry_agg_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, fill = sector_agg)) +
  geom_col() +
  facet_grid(rows = vars(branch), cols = vars(scenario), scales = "free_y") + 
  scale_fill_manual(values = pal_industry_sectors, name = "sector") + 
  labs(x = "", y = "Mt CO2", title = "Global industry CO2 emissions (no bio) by sector") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_emissions_no_bio_industry_sector_agg_scenario_branch_comp_global.png"), width = 30, height = 12)

Industry non-CO2 emissions

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg_total %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2e, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "Mt CO2e", title = "Total industry GHG emissions") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_total_scenario_branch_comp.png"), width = 30, height = 20)

ggplot(nonCO2_emis_sector_total_CO2eq_industry_agg_total_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2e, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "Mt CO2e", title = "Global total industry GHG emissions") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/GHG_emissions_industry_total_scenario_branch_comp_global.png"), width = 15, height = 8)

Industry water withdrawals

ggplot(industry_water_withdrawals_total_no_H2 %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~region, nrow = 4, scales = "free_y") + 
  labs(x = "", y = "km^3", title = "Total industry water withdrawals (excluding hydrogen sectors)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        axis.text.x = element_text(angle=55, vjust = 1, hjust = 1),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_total_scenario_branch_comp.png"), width = 30, height = 20)

ggplot(industry_water_withdrawals_total_no_H2_global %>% filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = value, linetype = branch, color = scenario)) +
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "km^3", title = "Global total industry water withdrawals (excluding hydrogen sectors)") + 
  theme_bw() + 
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/industry_water_withdrawals_total_scenario_branch_comp_global.png"), width = 15, height = 8)

Total carbon emissions

ggplot(CO2_emis_global_w_LUC_total %>%  filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, linetype = branch, color = scenario)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
   scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "Mt CO2", title = "Global total net CO2 emissions") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_net_scenario_branch_comp_global.png"), width = 15, height = 8)

ggplot(CO2_emis_global %>%  filter(year >= 2005 & !grepl("USA", scenario)),
       aes(x = year, y = Mt_CO2, linetype = branch, color = scenario)) + 
  geom_line(size = 1) + 
  scale_linetype_manual(values = c("dashed", "solid")) + 
  scale_color_manual(values = pal_scenarios) + 
  facet_wrap(~policy_type, ncol = 2) + 
  labs(x = "", y = "Mt CO2", title = "Global total FFI CO2 emissions") +
  theme_bw() +
  theme(text=element_text(size=17),
        strip.background = element_blank())

ggsave(paste0(fig_dir, "/CO2_FFI_scenario_branch_comp_global.png"), width = 15, height = 8)